#ifndef LUA_NUMERIC_FT_HPP
#define LUA_NUMERIC_FT_HPP
#include <omp.h>
#include <luaNumeric/config.h>
#include <vector>
#include <string>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <unsupported/Eigen/FFT>
#include <boost/bind.hpp>
#include <algorithm>
#include <numeric>
#include <complex>
#include <appendix/phy_constants.h>
#include <log/log.h>
#include <luaNumeric/vector.hpp>
namespace FT {

	//xx 是位置
	//wx是每一个位置上的场
	//kx是k空间波矢量
	//dx 实空间间隔
	template< class VecX, class VecWx >
	std::complex<double> fx2kx( const VecX& xx, const VecWx& vv, double kx, double dx )
	{
		if( (xx.size() != vv.size() ) )
			throw sol::error("size do not match in fx2kx function");
		std::complex<double> res(0);
		for(int ii=0; ii < xx.size(); ii ++)
		{
			res += dx * vv[ii] * std::exp( ci * kx * xx[ii] );
		}
		return res;
	}

	template< class VecX, class VecWx >
	void fx2kx( const VecX& xx, const VecWx& vv, std::vector<double> kx, double dx,  std::vector< std::complex<double> >& res )
	{
		if( (xx.size() != vv.size() ) ||
			(kx.size() != res.size()) )
			throw sol::error("size do not match in fx2kx function");

		//std::complex<double> res(0);
		#pragma omp parallel for
		for(int nkx=0; nkx < res.size(); nkx++){
			res[nkx] = 0;
			for(int ii=0; ii < xx.size(); ii ++){
				res[nkx] += dx * vv[ii] * std::exp( ci * kx[nkx] * xx[ii] );
			}
			//logs::info_log( "%g,%g",res[nkx].real(),res[nkx].imag());
		}
		//return res;
	}

	const double sqrt2pi = std::sqrt( 2.0 * pi );
	template< class VecX, class VecWx >
	std::complex<double> ft1D( const VecX& xx, const VecWx& vv, double kx, double dx )
	{
		if( (xx.size() != vv.size() ) )
			throw sol::error("size do not match in fx2kx function");
		std::complex<double> res(0);
		for(int ii=0; ii < xx.size(); ii ++)
		{
			res += dx * vv[ii] * std::exp( -ci * kx * xx[ii] ) / sqrt2pi;
		}
		return res;
	}

	template< class VecX, class VecWx >
	void ft1D( const VecX& tt, const VecWx& vt,
			std::vector<double> ww, double dt,  std::vector< std::complex<double> >& res )
	{
		if( (tt.size() != vt.size() ) ||
			(ww.size() != res.size()) )
			throw sol::error("size do not match in fx2kx function");
#pragma omp parallel for
		for(int nww=0; nww < res.size(); nww++){
			res[nww] = 0;
			for(int ntt=0; ntt < tt.size(); ntt ++){
				res[nww] += vt[ntt] * std::exp( -ci * ww[nww] * tt[ntt] );
			}
			res[nww] *= dt/sqrt2pi;
		}
	}

	template< class VecX, class VecWx >
	void ft1D( const VecX& tt, const VecWx& vt,
			   std::vector<double> ww, std::vector< std::complex<double> >& res )
	{
		if( (tt.size() != vt.size() ) ||
				(ww.size() != res.size()) ||
				(tt.size() < 2 ) )
			throw sol::error("size do not match in ft1D function");
		//double dt = tt[2]-tt[1];
#pragma omp parallel for
		for(int nww=0; nww < res.size(); nww++){
			res[nww] = 0;
			for(int ntt=1; ntt+1 < tt.size(); ntt ++){
				res[nww] += (tt[ntt+1]-tt[ntt-1])/2.0 * vt[ntt] * std::exp( -ci * ww[nww] * tt[ntt] );
			}
			res[nww] += (tt[2]-tt[1]) * vt[0] * std::exp( -ci * ww[nww] * tt[0] );
			res[nww] += (tt[tt.size()-1]-tt[tt.size()-2]) * vt[tt.size()-1]
														  * std::exp( -ci * ww[nww] * tt[tt.size()-1] );
			res[nww] /= sqrt2pi;
		}
	}

	template< class VecX, class VecWx >
	std::complex<double> ift1D( const VecX& xx, const VecWx& vv, double kx, double dx )
	{
		if( (xx.size() != vv.size() ) )
			throw sol::error("size do not match in ft1D function");
		std::complex<double> res(0);
		for(int ii=0; ii < xx.size(); ii ++)
		{
			res += dx * vv[ii] * std::exp( ci * kx * xx[ii] ) / sqrt2pi;
		}
		return res;
	}

	//频域变换到时域
	//ww: omega
	//vw：value of omega
	//tt： time
	//dw: delta omega
	template< class VecX, class VecWx >
	void ift1D( const VecX& ww, const VecWx& vw,
			std::vector<double> tt, double dw,  std::vector< std::complex<double> >& vt )
	{
		if( (ww.size() != vw.size() ) ||
			(tt.size() != vt.size()) )
			throw sol::error("size do not match in ift1D function");

#pragma omp parallel for
		for(int ntt=0; ntt < vt.size(); ntt++){
			vt[ntt] = 0;
			for(int nww=0; nww < ww.size(); nww ++){
				vt[ntt] += vw[nww] * std::exp( ci * tt[ntt] * ww[nww] );
			}
			vt[ntt] *= dw/sqrt2pi;
		}
	}

	template< class VecX, class VecWx >
	void ift1D( const VecX& ww, const VecWx& vw,
			    std::vector<double> tt, std::vector< std::complex<double> >& vtt )
	{
		if( (ww.size() != vw.size() ) ||
			(tt.size() != vtt.size())  ||
			(ww.size() < 2)
			)
			throw sol::error("size do not match in ift1D function");

#pragma omp parallel for
		for(int ntt=0; ntt < vtt.size(); ntt++){
			auto&& t = tt[ntt];
			auto& vt =  vtt[ntt];
			vt = 0;
			for(int nww=0; nww+1 < ww.size(); nww ++){
				auto wmid = (ww[nww+1]+ww[nww])/2.0;
				auto vwmid = (vw[nww+1]+vw[nww])/2.0;
				auto dw = ww[nww+1]-ww[nww];
				vt += dw * vwmid * std::exp( ci * t * wmid );
			}
			vt /= sqrt2pi;
		}
	}

	template< class VecT, class VecVT, class VecW, class VecVW >
	void conv1D( const VecT& tt, const VecVT& vt, double dt,  //输入的时间序列
			     const VecW& ww, const VecVW& vw, double dw,  //频域响应
				 const VecT& ttRes, std::vector< std::complex<double> >& vtRes)
	{
		if( (ww.size() != vw.size() ) ||
			(tt.size() != vt.size())  ||
			(ttRes.size() != vtRes.size()) )
			throw sol::error("size do not match in conv1D function");
		VecVW vwSrc(ww.size());
		ft1D(tt,vt,ww,dt,vwSrc);

#pragma omp parallel for
		for(int nw=0; nw < vwSrc.size(); nw++)
		{
			vwSrc[nw] *= vw[nw];
		}

		ift1D(ww,vwSrc,ttRes,dw,vtRes);
	}

	template< class VecT, class VecVT, class VecW, class VecVW >
	VecVW conv1D( const VecT& tt, const VecVT& vt,  double dt, //输入的时间序列
			      const VecW& ww, const VecVW& vw,             //频域响应
				  const VecW& ww1,double dw,                   //傅里叶变换区间
				  const VecT& ttRes, std::vector< std::complex<double> >& vtRes)
	{
		if( ww.size() != vw.size()  ||
			tt.size() != vt.size()  ||
			ttRes.size() != vtRes.size()||
			ww1.size() < 2 )
			throw sol::error("size do not match in conv1D function");
		//插值
		VecVW vwSrc1(ww1.size());
		ft1D(tt,vt,ww1,dt,vwSrc1);
		//ft1D(tt,vt,ww1,tt[2]-tt[1],vwSrc1);
		auto vw1 = Vector::interpolation_complex(ww,vw,ww1);
#pragma omp parallel for
		for(int nw=0; nw < vwSrc1.size(); nw++)
		{
			vwSrc1[nw] *= vw1[nw];
		}
		ift1D(ww1,vwSrc1,ttRes,dw,vtRes);
		return std::move(vwSrc1);
	}


	template< class VecT, class VecVT, class VecW, class VecVW >
	std::tuple<VecW,VecVW> conv1D( const VecT& tt, const VecVT& vt,  			//输入的时间序列
			     	 	 	 	   const VecW& ww, const VecVW& vw,             //频域响应
								   double wstart,double wend, int nw,           //傅里叶变换区间
								   const VecT& ttRes, std::vector< std::complex<double> >& vtRes)
	{
		if( (ww.size() != vw.size() ) ||
			(tt.size() != vt.size())  ||
			nw < 2 ||
			(ttRes.size() != vtRes.size()) )
			throw sol::error("size do not match in conv1D function");
		//插值
		VecVW vwSrc1(nw);
		VecW ww1 = Vector::linear_space(wstart, wend, nw);
		ft1D(tt,vt,ww1,vwSrc1);
		auto vw1 = Vector::interpolation_complex(ww,vw,ww1);
#pragma omp parallel for
		for(int nw=0; nw < vwSrc1.size(); nw++)
		{
			vwSrc1[nw] *= vw1[nw];
		}
		ift1D(ww1,vwSrc1,ttRes,ww1[2]-ww1[1],vtRes);
		return std::make_tuple(std::move(ww1), std::move(vwSrc1));
	}



	//xx yy是位置
	//vv是每一个位置上的场
	//kx ky 是k空间波矢量
	//ds 实空间面积
	template< class VecX, class VecWx >
	std::complex<double> ft2D( const VecX& xx, const VecX& yy, const VecWx& vv, double kx, double ky , double ds )
	{
		if( (xx.size() != yy.size() || xx.size() != vv.size() ) )
			throw sol::error("size do not match in ft2D function");
		std::complex<double> res(0);
		for(int ii=0; ii < xx.size(); ii ++)
		{
			res += ds * vv[ii] * std::exp( ci * kx * xx[ii] ) * std::exp( ci * ky * yy[ii] );
		}
		return res;
	}

	template< class VecX, class VecWx >
	void ft2D( const VecX& xx, const VecX& yy, const VecWx& vv,
			std::vector<double> kx,
			std::vector<double> ky,
			double ds,  std::vector< std::complex<double> >& res )
	{
		//std::complex<double> res(0);
		if( (xx.size() != yy.size() || xx.size() != vv.size() ) ||
			(kx.size() != res.size() || ky.size() != res.size() ) )
			throw sol::error("size do not match in ft2D function");
		#pragma omp parallel for
		for(int nkx=0; nkx < res.size(); nkx++){
			res[nkx] = 0;
			for(int ii=0; ii < xx.size(); ii ++){
				res[nkx] += ds * vv[ii] * std::exp( ci * kx[nkx] * xx[ii] ) * std::exp( ci * ky[nkx] * yy[ii] );
			}
		}
		//return res;
	}

	template< class VecX, class VecWx >
	void ft2D( const VecX& xx, const VecX& yy, const VecWx& vv,
			std::vector<double> kx,
			std::vector<double> ky,
			const VecX& ds,  std::vector< std::complex<double> >& res )
	{
		//std::complex<double> res(0);
		if( (xx.size() != yy.size() || xx.size() != vv.size() || xx.size() != ds.size() )
			||(kx.size() != res.size() || ky.size() != res.size() ) )
			throw sol::error("size do not match in ft2D function");
		#pragma omp parallel for
		for(int nkx=0; nkx < res.size(); nkx++){
			res[nkx] = 0;
			for(int ii=0; ii < xx.size(); ii ++){
				res[nkx] += ds[ii] * vv[ii] * std::exp( ci * kx[nkx] * xx[ii] ) * std::exp( ci * ky[nkx] * yy[ii] );
			}
		}
		//return res;
	}



	void efft(Eigen::MatrixXd& in, Eigen::MatrixXcd& out);

};

#endif
